AD- A 188  473 


VT'.T't'  A.Ttyx.T  ’.T  'TTT? 


FINAL  REPORT 


T"  ’JL"  1  -L"  T '  V  -J  %  'VT1IW1  ^ 


($1 


OIK  F/L£  copy 


•*.V 

v.*;, 

;.*vi 

■.-vl 

t  V*  ^ 


AN  EFFICIENT  NUMERICAL  ALGORITHM  FOR 
SOLVING  SCATTERING  AND  INVERSE  SCATTERING 
PROBLEMS  OF  ELECTROMAGNETIC  WAVES 


Prepared  by  Y.  M,  Chen 
NUMERICAL  COMPUTATION  CORP. 

57  Quaker  Path 

Stony  Brook,  New  York  11790-1309 
Phone:  516-751-9518 

September  21,  1987 


APPROVED  FOR  PUBLIC  RELEASE 
DISTRIBUTION  UNLIMITED 


i mi 


l  >: 
, 


DTIC 

D 

SBIR  PHASE  I  Contract  No.  N00014-86-C-0109 
Office  of  Naval  Research 
Department  of  the  Navy 
800  N.  Quincy  Street 
Arlington,  VA  22217-5000 

SCIENTIFIC  MONITOR: 

Defense  Advanced  Research  Projects  Agency  -  TTO 
1400  Wilson  Blvd. 

Arlington,  VA  22209 


BLECTE 
0EC3  OK*r 


r 


87  12  il  108 


mm 


REPORT  DOCUMENTATION  PAGE 


la.  REPORT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 


2a.  SECURITY  CLASSIFICATION  AUTHORITY 
Directorate  for 


2b.  DECLASSIFICATION  /  DOWNGRADING  SCHEDULE 


Ul'il-i.'IMviKJ  HI  /»! 


4.  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 
FINAL  REPORT 


6a.  NAME  OF  PERFORMING  ORGANIZATION 
NUMERICAL  COMPUTATION  CORP. 


6c  ADDRESS  (City,  State,  and  ZIP  Code) 

57  Quaker  Path 

Stony  Brook,  N.Y.  11790-1309 


8a.  NAME  OF  FUNDING /SPONSORING 
ORGANIZATION 

Office  of  Naval  Research 


8c.  ADDRESS  (City,  State,  and  ZIP  Code) 
800  N.  Quincy  Street 
Arlington,  VA  22217-5000 


lb  RESTRICTIVE  MARKINGS 


3.  DISTRIBUTION  /AVAILABILITY  OF  REPORT 
APPROVED  FOR  PUBLIC  RELEASE 

DISTRIBUTION  UNLIMITED 


S.  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 
FINAL  REPORT 


7a.  NAME  OF  MONITORING  ORGANIZATION 
DARPA  -  TTO 


7b  ADDRESS  (City,  State,  and  ZIP  Code) 
1400  Wilson  Blvd. 
Arlington,  VA  22209-2308 


8b.  OFFICE  SYMBOL  9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 
(If  applicable) 

N00014-86-C-0109 


6b.  OFFICE  SYMBOL 
(If  applicable) 


10.  SOURCE  OF  FUNDING  NUM8ERS 


PROGRAM 
ELEMENT  NO. 


PROJECT 

NO. 


WORK  UNIT 
ACCESSION  r. 


11.  TITLE  (Include  Security  Classification) 

An  Efficient  Numerical  Algorithm  for  Solving  Scattering  and  Inverse  Scattering 
Problems  of  Electromagnetic  Waves  -  Unclassified 


12.  PERSONAL  AUTHOR(S) 


13a.  TYPE  OF  REPORT 
Final 


Chen,  Yung  M. 


13b.  TIME  COVERED  1 14.  DATE  OF  REPORT  (Year.  Month,  Day)  IlS.  PAGE  COUNT 

from  85- 12- 15to  87-8-31 1  87-9-21  I  44 


17. _  COSATI  COOES _ I  18.  SUBJECT  TERMS  (Continue  on  reverse  if  necessary  and  identify  by  block  number ) 

FIELD  |  GROUP  |  SUB-GROUP  |  „ _ .  _ _ 


Numerical,  algorithm,  scattering,  electromagnetic, 
wave 


19  ABSTRACT  (Continue  on  reverse  if  necessary  and  identify  by  block  number)  *£he  development  of  an  efficient  numerical 
algorithm  of  determining  the  unknown  material  composition  and  shape  of  an  arbitrary  target 
(from  the  measured  electromagnetic  waves  in  the  far  field  region  will  enhance  the  capability 
of  the  defense  radar  system  to  defeat  known  evasive  schemes.  The  first  step  in  this  researc 
effort  is  the  development  of  an  efficient  and  versatile  numerical  algorithm  for  calculating 
the  scattered  electromagnetic  waves/radar  cross  section  by  a  target  with  known  complex 
geometry  and  material  property.  Hence  the  pupose  of  this  Phase  I  research  is  to  develop  an 
efficient  numerical  algorithm  for  solving  twordimensional  scattering  problems.  This  is 
achieved  by  using  a  special  finite  difference  method  based  upon  a  natural  spatial  discreti¬ 
zation  of  the  integral  form  of  Maxwell's  equations  on  a  non-orthogonal  grid-system  and  the 
leap-frog  finite  differencing  in  the  time  domain.  It  has  the  advantages  of  being  (a)  more 
efficient  than  any  other  known  numerical  methods,  (b)  highly  accurate  due  to  the  body-fitted 
grid  system,  and  (c)  the  easiest  numerical  method  to  implement  boundary  conditions.  The 
capability  and  feasibility  of  this  two-dimensional  computer  code  are  tested  by  performing - 


20.  DISTRIBUTION  /AVAILABILITY  OF  ABSTRACT 
(3 UNCLASSIFIED/UNLIMITED  El  SAME  AS  RPT. 


22*.  NAME  OF  RESPONSIBLE  INDIVIDUAL 


121.  ABSTRACT  SECURITY  CLASSIFICATION 

Unclassified 


22b.  TELEPHONE  (Include  Ana  Code) 


00  FORM  1473,84  MAR 


83  APR  edition  may  bt  used  until  exhausted. 
All  other  editions  are  obsolete. 


v-  Block  19: 

V 

numerical  simulations  on  few  realistic  examples,  e.g. ,  cylindrical  objects  with 
cross  sections  of  a  metallic  jet  and  a  composite  airfoil.  In  these  processes, 
the  radar  cross  sections  as  functions  of  both  the  incident  angle  and  the  scattered 
angle  are  calculated  and  they  seem  to  be  quite  good.  Hence  this  Phase  I  research 
is  completed  successfully.  <•  .  ps  f.'  ^  V  >i  -  c 

>  !  (  '  "  '  *<-  A  u  <->  >  ■  c  ^CA  a-  ^ 


SBIR  PHASE  I  CONTRACT  NO.  N00014-86-C-0109 

OFFICE  OF  NAVAL  RESEARCH 

DEPARTMENT  OF  THE  NAVY 

800  N.  Quincy  Street 

Arlington,  Virginia  22217-5000 

TITLE:  AN  EFFICIENT  NUMERICAL  ALGORITHM  FOR  SOLVING  SCATTERING 
AND  INVERSE  SCATTERING  PROBLEMS  OF  ELECTROMAGNETIC  WAVES 

FINAL  REPORT:  Prepared  by  Y.  M.  Chen 

NUMERICAL  COMPUTATION  CORP. 

57  Quaker  Path  Ph.  516-751-9518 

Stony  Brook,  N.Y.  11790-1309 


DATE:  September  21,  1987 


TABLE  OF  CONTENTS 

SUMMARY 

DEGREE  TO  WHICH  PHASE  I  OBJECTIVES  HAVE  BEEN  MET 
ANTICIPATED  BENEFITS 
INTRODUCTION 

FINITE  DIFFERENCE  METHOD 
NUMERICAL  SIMULATION 
TECHNICAL  DISCUSSION 
REFERENCES 

APPENDIX  -  COMPUTER  CODE  ,  .  ^ 

Justification. 


Accession  Tor 


WXS  OtAAX 

MIC  fAl 


□ 

□ 


Distribution/ 


Availability  Cedes 


Dist 

/\ 


* 


{Avail  and/or 
Spsolal 


i 

ii 

ii 

1 

4 

8 

11 

19 

20 


SUMMARY 


The  development  of  an  efficient  numerical  algorithm  capable  of  determining 
the  unknown  material  composition  and  shape  of  an  arbitrary  target  from  the 
measured  electromagnetic  waves  in  the  far  field  region  will  enhance  the  capabi¬ 
lity  of  the  defense  radar  system  to  defeat  known  evasive  schemes.  The  first 
step  in  this  research  effort  is  the  development  of  an  efficient  and  versatile 
numerical  algorithm  for  calculating  the  scattered  electromagnetic  waves/radar 
cross  section  by  a  target  with  known  complex  geometry  and  material  property. 
Hence  the  purpose  of  this  Phase  I  research  is  to  develop  an  efficient  numerical 
algorithm  for  solving  two-dimensional  scattering  problems.  This  is  achieved 
by  using  a  special  finite  difference  method  based  upon  a  natural  spatial 
discretization  of  the  integral  form  of  Maxwell's  equations  on  a  non-orthogonal 
grid-system  and  the  leap-frog  finite  differencing  in  the  time  domain.  It  has 
the  advantages  of  being  (a)  more  efficient  than  any  other  known  numerical 
methods,  (b)  highly  accurate  due  to  the  body-fitted  grid  system,  and  (c)  the 
easiest  numerical  method  to  implement  boundary  conditions.  The  capability 
and  feasibility  of  this  two-dimensional  computer  code  are  tested  by  performing 
numerical  simulations  on  few  realistic  examples,  e.g. ,  cylindrical  objects 
with  cross  sections  of  a  metallic  jet  and  a  composite  airfoil.  In  these 
processes,  the  radar  cross  sections  as  functions  of  both  the  incident  angle 
and  the  scattered  angle  are  calculated  and  they  seem  to  be  quite  good. 

Hence  this  Phase  I  research  is  completed  successfully. 
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DEGREE  TO  WHICH  PHASE  I  OBJECTIVES  HAVE  BEEN  MET 

The  purpose  of  the  Phase  I  research  is  to  develop  an  efficient  numerical 
algorithm  for  solving  two-dimensional  scattering  of  electromagnetic  waves  by 
a  target  with  complex  geometry  and  material  property.  Now,  we  are  happy  to 
announce  that  the  objectives  of  the  Phase  I  research  have  been  achieved 
completely  and  successfully.  The  complete  2-D  computer  code  is  given  in  the 
Appendix  in  this  final  report. 


ANTICIPATED  BENEFITS 

As  it  stands  by  itself,  the  usefulness  of  the  2-D  computer  code  developed 
here  in  Phase  I  research  is  rather  limited,  because  in  the  real  life  there  is 
no  truly  two-dimensional  target.  However,  it  has  demonstrated  beyond  any 

doubt  that  the  above  developed  numerical  algorithm  can  be  extended  to  compute 
efficiently  and  accurately  the  scattered  E-M  waves  or  radar  cross  section 
for  the  three-dimensional  target  with  complex  geometry  and  material  property. 
The  3-D  code  will  be  an  important  subroutine  of  the  Generalized  Pulse- 
Spectrum  Technique  for  solving  the  3-D  inverse  scattering  problems;  moreover, 
it  will  be  a  very  useful  tool  in  advanced  vulnerability  analysis  for  any 
postulated  target. 
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INTRODUCTION 


The  numerical  methods  for  solving  the  electromagnetic  waves  scattering 
problems  have  a  recent  history  of  only  two  decades.  Although  the  true  and 
original  radar  scattering  problem  is  a  problem  of  solving  the  initial-boundary 
value  problem  of  Maxwell's  equations  in  the  space-time  domain,  in  earlier  days 
most  numerical  methods  are  based  upon  first  Fourier  transforming  the  original 
problem  in  the  space-time  domain  into  a  corresponding  problem  in  the  space- 
frequency  domain  and  then  solving  it  numerically  for  several  chosen  frequencies. 
In  this  way,  some  useful  but  incomplete  scattering  information  can  be  extracted 
from  these  results  without  performing  inverse  Fourier  transformation.  In 
particular,  the  antique  definition  of  the  radar  cross  section  was  originally 
defined  for  a  single-frequency  source;  its  proportionality  to  the  square  of 
the  electric  field  and  the  invalidity  of  the  principle  of  superposition  for 
this  nonlinear  situation  will  introduce  further  error  into  the  calculation  of 
radar  cross  section  by  this  approach  as  compared  to  the  true  situation. 

Moreover,  all  numerical  methods  for  solving  the  E-M  scattering  problem  in  the 
space-frequency  domain  can  be  reduced  to  the  problem  of  solving  a  very  large 
matrix  equation  (in  particular  the  three-dimensional  scattering  problems) 
which  requires  an  extra  ordinarily  large  amount  of  CPU  time  and  memory  storage; 
thus  it  renders  these  methods  inefficient. 

An  efficient  finite  difference  method  for  solving  E-M  scattering  problems 
in  the  space-time  domain  was  first  introduced  by  Yee  [l]  for  the  two-dimensional 
case  and  later  applied  to  the  three-dimensional  case  by  Taflove  and  Brodwin  [2] , 
Holland  [3] ,  and  Kunz  and  Lee  [4] .  This  finite  difference  method  requires 
a  uniformly  rectangular/cubical  grid  system  and  it  is  the  most  efficient  method 
(see  the  review  by  Chen  [5]).  Unfortunately,  for  scatterers  with  curved 
boundaries,  one  needs  extra  ordinarily  large  amount  of  uniform  rectangular/ 
cubical  grid  zones  to  approximate  the  curved  boundaries  and  minimize  the 
undesirable  "staircasing  phenomenon",  and  thus  it  renders  this  numerical  method 
inefficient  in  general. 

Later,  Mei,  Cangellaris,  Angelakos  and  Lin  [6] ,  [7]  have  presented  the 
"Point-matched  Time  Domain  Finite  Element  Method",  a  combination  of  the 
essential  features  of  the  standard  Yee's  finite  difference  method  and  the 
finite  element  method.  Its  efficiency  is  an  improvement  over  the  standard 
Yee's  finite  difference  method  due  to  its  ability  to  compute  on  a  body-fitted 


non-orthogonal  grid  system,  but  it  is  still  not  efficient  enough  due  to  the 
additional  work  in  implementing  any  boundary  conditions,  e.g.,  extra  inter¬ 
polation  and  extrapolation  are  needed  at  the  boundary  grid  zones. 

Recently,  Yee  C8]  has  made  a  dramatic  improvement  of  his  method  by  applying 
his  finite  difference  discretization  in  the  most  natural  way  to  the  integral 
form  of  Maxwell's  equations  on  a  general  non-orthogonal  grid  system  which 
makes  the  implementation  of  boundary  conditions  extremely  easy.  One  can 
show  that  it  is  the  most  efficient  (for  the  same  accuracy)  method  available 
by  performing  the  standard  computational  complexity  analysis,  i.e.,  to  count 
the  total  floating  point  arithmetic  operations  needed  in  a  typical  calculation 
[5].  This  algorithm  can  be  vectorized  and  parallelized  with  great  ease  and 
hence  it  will  be  ideal  for  the  vector  and  multi-processor  computers. 

Here  in  the  Phase  I  research,  Yee's 
improved  method  is  generalized  for 
solving  the  two-dimensional  scattering 
problems  of  E-M  waves  by  targets  with 
complex  geometry  and  material  property. 

First,  the  whole  space  domain  fl  is 
divided  into  three  connected  but  non¬ 
overlapping  sub-domains,  the  interior 
region  representing  the  target 
and  possessing  a  non-orthogonal 
cylindrical  grid  system  centered  in 
itself,  the  intermediate  region  ^ 
representing  the  free  space  just 
outside  of  the  target  and  possessing 
the  same  grid  system,  and  the  exterior 
region  representing  the  far- field 
space  but  truncated  at  a  large  distance  away 

from  the  target  and  possessing  the  standard  Fig.  1 

orthogonal  cylindrical  grid  system  (Fig.  1). 

Mathematically,  an  initial-boundary  value  problem  of  the  integral  form  of 
Maxwell's  equations  must  be  solved  numerically, 
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f  E.  •  d£  =  -  £f  jj3H  /3t  •  ds, 

“  x  e  fl.,  (1) 

<j>  Hx  •  d£  =  #  (0E1  4-  eSEj/St)  •  ds, 

<f>  E,  •  d£  =  -  #  Mn3H,/3t  •  ds, 

L  0  L  xtSl,  (2) 

^  H2*di  =  ^  e0aI2/3t  *  d®* 

^  E  •  d£  =  -  3H  /3t  •  ds, 

J  °  J  x  e  0  ,  (3) 

^  I3  •  d£  =  #  ( J  +  e03E3/3t)  •  ds, 

with  boundary  conditions  (assuming  no  surface  charges  and  currents), 

x  e  3n12»  (4) 


(5) 

where  n  is  the  unit  outer  normal  vector  at  the  interfaces,  3fl12  is  the 
interface  between  ^  and  is  the  outer  boundary  of  Jlj,  J(x)  is  the 

source  distribution,  €q  and  are  the  free  space  permittivity  and  permeabi¬ 
lity  respectively,  £(x)  is  the  3  x  3  real  positive  symmetric  permittivity 
matrix  of  the  target,  y(x)  is  the  3  x  3  real  positive  symmetric  permeability 
patrix  of  the  target,  and  £(x)  is  the  3  *  3  real  positive  symmetric  conduct¬ 
ivity  matrix  of  the  target. 

Since  there  exists  no  realistic  two-dimensional  E-M  scattering  problem, 

for  the  Phase  I  research  the  scattering  of  normal  incident  TEM  electromagnetic 

wave  by  a  cylindrical  target  with  its  axis  along  i  ,  e.g. ,  !  ■  E  i  +E  i 

—z  x— x  y—y 

and  H  ■  H  i  ,  where  i  is  the  unit  vector  in  the  a-direction. 

—  z—z  —a  — 


n  x  E2  =  n  x  e  ,  n  x  H2  =  n  x  , 

and  the  asymptotic  terminating  condition, 

n  x  E3  =  (^q/ e0)**(jl  *  H.3> ,  x  E  3«3  , 
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FINITE  DIFFERENCE  METHOD 

To  discretize  (l)-(3),  the  rectangle  rule  is  used  to  approximate  both 
the  line  and  area  integrals,  the  values  of  12  and  H  fields  are  calculated  on 
two  different  but  staggered  grid  systems,  and  the  leap-frog  finite  difference 
scheme  is  used  to  approximate  the  first  order  derivative  in  time. 

Let  each  grid  point  of  the  basic  non-orthogonal  cylindrical  grid  system 

be  denoted  by  (r  .,  0  )  s  (i,j),  where  "i"  and  "j"  denote  the  i-th  closed 
i»J  J 

cylindrical  grid  line  and  the  j-th  radial  grid  line  respectively.  Let  the 
center  of  the  quadrilateral  defined  by  (i,j),  (i+l,j),  (i,j+l)  and  (i+l,j+l) 
be  denoted  by  (i+*S,j+*s)  ■  ri+i ,  j+ ri,  j+i+ ri+l ,  j+I>  •  ^  6j+l+ 6  J  ’ 

where  all  these  centers  form  another  non-orthogonal  cylindrical  grid  system 
staggered  on  the  basic  grid  system.  Let  the  j2  fields  be  evaluated  at  the 
mid-points  of  the  four  edges  of  the  quadrilaterals  (Fig.  2)  and  at  the  integer 

time  increments  t  =  nAt,  n  =  1,2,3, . ,;  let  the  H  fields  be  evaluated  at 

the  centers  of  the  quadrilaterals  (Fig.  2)  and  at  the  half  time  increments 
t  =  (n+*s)At,  n  *«  0,1, 2, 3 .  For  simplicity,  let 


€  »  0  e_ 


fy  0  0 
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u  =  0  y  0 
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source  located  at  the  grid  point 
(a, 6).  Let  i  *  1,2,3,... ,1 
and  j  =  1,2,3,. ..,J.  v/ij. 


a  **  a(a  scalar)  ,  and  J  be  a  point 
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Fig.  2 
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The  discretization  of  (l)-(3)  in  are, 

2At 


=  Hft.lWl.n-H)  -  ^  j+y  (S^-e^r^  .  tj “  - 

(E0s,j,n)r  +  E(l,j+4,n)6£(l,j+%)  -  E(h,  j+1  .n)^  ^  > , 

{Ata(Js,  j)6£2(%,  j)  +  e,  (^.j)  }E(%,  j  ,n+l)  (7) 


-  Ato(lSij){r.  (l-(0.+1-e  )^/S)  -  r,  .  (l-(6  -0  .)*/8)}6£(%,j)J  1  £  E^.j.n) 

'i.J+'S  J+1  J  'StJ~'S  J  j 

+  e,  (*5,  jJEte.j  ,n)  +  At6£*(JS,j)  (H(%,  j+^.n+y  -  H(*s,  j-^.n+^J )  ,  j  =  1,2,3,... ,J, 

n  =  0, 1,2,3, . . 

8(1+^, ]+h,n+h)  -  H(i+J<j, j+%,n-is)  (8) 


At 


Pz(i+JS,j+Js)5A(i+%,j+Js){E(l+l2,:i  ,n)  (ri+l,j“ri,j) 

+  E(i+l,j+%,n)6£(i+l,j+%)  -  E(i+%, j+1 ,n) (r±+1  j+1) 

-  E(i,  j+*s,n)6£(i,  j+^s)  } , 


{(Ata(i,j+Js)+ea(i,j4^))6£(i,j445  }E(i,  j+Js,n+l)  (9) 

-  {Ata(i,j+Js)(ri  J+1-ri}j)(l-<ej+1-ej)2/8)  +  eg(i,j+is) }E(i-^,j  ,n+l) 

=  6i(i,j+Js)ea(i,j+ii)E(i,j+J5,n)  -  eg(i,  j+Js)E(i+Js,  j  ,n) 

-  —■ -  (H(i+Hj,j+Js,n44s)  -HCi-Js.j+^.n+^s))  , 

i+,S,j+iS  i-JS,j+JS 

-{4to<i4l,,J)[r1^(;J^(l-(ej+1-8j)2/8)-r1+1!l>jJj(l-(ej-0j.I)2/8)l«c(i-«,.J)}-  (10) 

6£(i+ij,  j)E(i,  j-Ms,n+l) 

+  {Ata(i+*s,  j)SS,2(i+*s,  j)  +  ex(i-M5,j)}E(i+S5,j,n+l) 

-  -  fi^d+Js.jJe^U+^.jJEa.j+is.n)  +  ex(i+Js,j)E(i+Js,j,n) 

+  At6£*(i+ls,j)(H(i+ls,j+5s,n+Js)  -  H(i+ii,  j-^n-^)  )  , 

j  *  1,2,3 . J,  n  ■  0,1, 2, 3 . 


The  discretization  of  (l)-(3)  in  are» 
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H(i+*s,  j+*s,n+*s)  -  H(i+*s,j+Ji,n-ls) 

_  At 

Jo  1 

-E(i+Ji,j+l,n)(ri+1  j+1-ri>j+1)  -  E(i,j+»i,n)6l(i,  j+ij)}, 
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{  (AtaQ  +  e0)<5£(i,  j+^)  }E(i, j+^,n+l) 

■  (rij+rrij)(1-(Vr0j)2/8)(At0o+eo)E(i+is’j’tl+1) 

=  <SS-(i, j+lj) EqECI, j+%,n)  -  eo(rijj+1-rijj)(l-(0j+1-0.)2/8)E(i+J5,j,n) 

_  - Wi-Hs.j-Us.n-^)  -  HU-h,i+h,n+h)), 

-  (4w0+e0){r1+lj >J+>s(1-(0j+1-ej,2/8)-ri+!i>J_l5(i-(e.-e._1)2/8)}«j(i+ii,j)E(i,j+|s,„+i) 
+  (Ato0+e())6<l2(i+Js,j)E(i+ls,j,n+l)  (13) 

+  e^U+^.JJEd+Js.J.n)  +  Ac6^*(i+ii,  j)  (H(i+is,  j+Jjtn+!s)  - H(i+!j, j-*s,n+*s) ) 
j  =  1,2,3,. ..,J,  n  =  0, 1,2,3, . 


The  discretization  of  ( 1) — ( 3)  in  are» 

H(i+*i,j+Js,n+ls)  =  H(i+Js,  j+!$,n-%) 

2At 


(14) 


~  7"(e - Z0-)('r - -r  )  (r - ,n)-E(i+Js,j+l,n))  (r  -r  ) 

j+1  i+1  iM  i+1  V  1+1  1 

+  (ri+1E(i+l,j44s,n)  -r.E(i,j+l5,n)X0:)+1-0j)}, 

E(i,j+is,n+l)  -  (Ata  +e  )_1{e  E(i,j+>s,n)  -  2At  - v  (H(i-Hj, i+%,n+h)-n(i-H,  j-Ms.n-Hj) ) 

v  i+1  i-l' 

~  3At  (  ~  Jx(a,6,n+l5)sinl5(ej+1+0_j)  (15) 

+  Jy(a,6,n+%)cosJ5(6j+1+e  J)}, 

eu^,+o  - 

-<$ i , a6  j+*s ,  0A 1 C" Jx (° » e  * n44s)  ( s inJs( 6 j + 1+0 j ) _s inJ* ( 0 j' +e j  _ i ) ) 

+Jy(a,0,n+»s)  (coslj(0j+1+8j)-cosJ5(ej+ej_1))]},  (16) 

j  =  1»2,3 . J,  n  “  0,1, 2, 3, . , 

where  6  A  (!+*,;)+»*)  -  W0J+1-0j){(r1+l#J-r1|J)r1#J+1  +  (*1+1 

6**(i,J+*s)  - 

s^d+h.j)  -  ^'i^.jWrV  +  ri+Js,j-Js^8j~0j-i^’ 

ri+Jj,j-Mj  “  *s(ri,j+ri,j+l+ri+l,j+ri+l,j+l)' 
ri+iS,j-!3  "  J‘(ri,j-l+ri,j+ri+l,j+ri+l,j)' 
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<$£(i,  j+*i) 
5Jl(i+Js,  j) 
ea(i»j+J5) 
Eg(i,j+%) 

eA  (i+is,  j) 


e?(i+l5j) 


ex(i,j+%) 

ey(i,j+»s) 

ex(i^a.j) 

ey(i+Js,j) 


Ur, 


'-r<  <>2  + 


i, j+1  i, j  u 

^ri+i5,j+lS  ri+JS,j-i5^  +  !iri+l5, j+1sri+5S, 9j-l^  J  ’ 
v(i»j+4)sin\(0  .+9  )  +  e„(i,j+5s)cos2is(0  .+0 .), 

x  j+i  j  y  j+i  j 


9J+1  -  t^.sin  e^.inWe^-Hy 
+  ey(i.J+y(r1J+|cos  0.+1  -  r^cos  6 .  JcosW^+S^ , 
ex(i+ls.l);r^  +  r2^  j^sin2!^  .+0 

_2ri+*s.3+%,ri+‘s.J-%sitlls<ej+l+ej)sIn,iCej+ej-l)} 

+  e^i+b.jHr^^cos^te^^)  + 

■2ri^o+>iri^,  J->sC03ls<9j+i+83)c0S’i(ej+ej-i) }  • 

£x(i-^,j)si"  VrWi,rt‘ln!i(Vl+V  "  ri+>l,j->iSinlS<6j+9j-l)) 

+  Eya+4.j)C0S  ej  (ri+>5,J+lsC°s!*<ej+l+0J>  -  ri+ls,j-liCOS'5<ej+0j-l))- 


>s(e  (i+^.j+Js)  +  e  (i-%,j+Js)), 

X  X 


5i(ey(i+Ss,j+ij)  +  EyCi-^J+is)), 


^(e  (i+*s,j+*s)  +  ej{(i-»45»j-!s)) , 


Js(£y(i+^s»  j+*i)  +  EyCi+%,  j-*s)) 


Theoretically,  the  boundary  conditions  (4)  must  be  imposed  at  the  inter¬ 
face  of  two  different  materials.  But  here,  there  is  no  need  to  impose  the 
boundary  conditions  explicitly  in  programming  this  numerical  algorithm, 
because  the  boundary  condition  for  the  tangential  component  of  12  is  satisfied 
automatically,  and  the  other  three  boundary  conditions  are  also  automatically 
satisfied  in  the  approximate  sense  if  the  differences  of  the  material  properties 
spread  linearly  across  a  complete  grid  zone  instead  of  just  across  the  inter¬ 
face.  In  this  way,  there  is  no  cumbersome  programming  instruction  at  the 
interface  to  slow  down  the  calculation  and  the  application  of  boundary 
conditions  is  replaced  by  the  process  of  assigning  the  material  parameters  into 
the  grid  zones.  In  particular,  the  programming  Instruction  is  extremely 
simple  if  the  material  interface  is  located  either  at  a  constant  i-line  or  at 
a  constant  j-line. 
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J 


Ac  the  exterior  boundary  (i  =  1) »  the  following  simple  but  effective 
discretization  of  the  asymptotic  terminating  of  non- re flee ting  condition  (5) 


is  used, 


E(I,j+is,n+l)  =  (M0/e0)'SH(I-^,j+i5,n+%), 


j  =  1,2,3,... ,J,  n  =  0, 1,2,3,, 


Finally,  to  achieve  a  stable  computation  for  this  explicit  finite 
difference  scheme,  one  must  impose  the  following  well  known  stability  condition. 

At  <  Min.  6J,  •  Min.  e  _(i+%,j+%)  •  Min.  u  „  (i,  j ) ,  (18) 

i,j  i, j ,x,y , z  ,y’  i,j,x,y,z  ,y’ 

where  52^  is  the  typical  dimension  of  the  i-jth  grid  zone. 


NUMERICAL  SIMULATION 


From  Eqs.  (6)-(16),  the  discrete  values  of  the  E-M  fields  in  the  space- 
time  domain  scattered  by  a  two-dimensional  target  with  complex  geometry  and 
material  property  can  be  calculated  as  accurate  as  one  desires.  However, 
the  main  interest  for  the  radar  scattering  application  is  not  the  detail 
description  of  the  E-M  fields  everywhere;  rather  it  is  a  certain  average 
quantity  to  characterize  the  scattering  and  absorbing  properties  of  a 
scatterer  in  the  incident  E-M  fields.  The  quantity  mostly  common  used  for 

this  purpose  is  the  radar  cross  section  which  initially  was  defined  for  the 
time-harmonic  E-M  fields  [9]  as  a  normalized  scattered  power  intensity  averaged 
over  a  cycle  or  a  normalized  scattered  energy  intensity  per  cycle, 

I  n  I  2 


lim  T 


UjJ2  ' 


where  T  *  2nr  for  two-dimensional  problems, 

2 

■  4irr  for  three-dimensional  problems, 

r  -  distance  between  the  transmitter  and  the  scatterer,  • 

i p  ■  angle  between  the  transmitter  and  the  zero  angle  axis, 

E°  ■  Incident  E-field  at  the  target, 

—in 

and  E  ■  scattered  far  field  E-field. 


Since  the  antique  definition  of  the  radar  cross  section  (19)  was  defined 
only  for  a  single-frequency  source,  errors  will  be  introduced  when  it  is 
applied  to  the  scattering  problems  with  general  time-varying  sources.  In 
order  to  be  able  to  deal  with  more  realistic  situations,  the  definition  of  the 
radar  cross  section  (19)  is  generalized  to  cases  with  C-W  pulse  as  source  as 
the  normalized  scattered  energy  intensity, 

/Ts+ATs|e  I2  dt 

ZW  =  lim  r  Ts  _  ,  (20) 

r-*»  t.  4AT .  0 

/  ln  ln|E°  |2  dt 

1  ™  I  _ln  I 

in 

where  AT.  is  the  duration  of  the  incident  E-M  pulse  and  AT  is  the  duration 
in  s 

of  the  scattered  E-M  pulse,  i.e.,  E.  5  0  for  t  <  T,  and  t  >  T.  +AT.  ,  and 

r  -=ln  in  in  in 

E  s  0  for  t  <  T  and  t  >  T  +AT  . 

— s  s  s  s 

It  is  clear  that  the  definition  (19)  is  a  special  case  of  the  definition  (20). 


Now,  there  are  two  approaches  to  calculate  the  radar  cross  section  of  a 
given  target.  One  approach  is  first  to  approximate  the  time-harmonic  source 
by  a  C-W  pulse  with  very  long  duration,  next  one  runs  the  above  developed 
computer  code  long  enough  for  all  of  the  transients  to  died  off,  and  finally 
the  definition  (19)  is  used  to  compute  the  radar  cross  section.  The  other 
approach  is  to  use  a  very  short  C-W  pulse  as  the  source,  next  the  computer 
code  is  run  long  enough  for  the  scattered  C-W  pulse  passing  through  the  loca¬ 
tion  of  the  receiver  completely,  and  finally  the  definition  (20)  is  used  to 
compute  the  radar  cross  section.  It  is  obvious  that  the  cost  for  the  first 
approach  to  calculate  the  radar  cross  section  is  prohibitively  high  and  the 
second  approach  is  the  preferred  one. 

To  achieve  our  goal,  a  simple  and  almost  costless  sub-routine  is  added  to 

the  computer  code  for  calculating  /  |E(x,T)| 2dl  in  the  center  of  everv  grid 

0 

zone  at  the  outer  boundary  dfty  i.e.  ,  ♦ 

l  |E(I-is,j+»s,m)rAt,  j  -  1,2,3,.  ...J,  n  ■  1,2,3, . . 

m*l 

are  calculated  in  the  code;  in  general,  it  is  a  monotonic  non-decreasing 
function  of  t. 


In  the  region  where  the  durations  of  E,  and  E  pulses  are  distinctively 
non-overlapping  (t|/  <  100°) ,  the  value  of  the  energy  integral  rises  to  a 
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T .  +AT .  2 

constant  value  /  in  xn  j  E^C >c , T )  |  dT  after  the  passing  through  of  the  the 

0  .Tg+ATg 

pulse  and  later  it  will  rise  again  to  another  constant  value  J  |E(3c,x)|zdT 

after  the  passing  through  of  the  E  pulse.  Then  the  value  of 

I  1 2  .  ~ S.  ,  .Ts+Ats,  1 2  /•^■in+^i 

jE^I  dl  can  be  obtained  simply  from  |I£|  dT  ~  Jq 

A  schematic  diagram  of  this  is  given  in  Fig.  3.  „  , 

T.j  +ATj  2  T-  +^Ti 

As  for  the  value  of  /  n  n|E?  |  dT,  it  can  be  obtained  as  /  Xn  n|E 


/ 


ts+ATS| 


I  E I 2dT. 


Lm 


-in' 


—in 1 


'dT 


at  the  target  location  with  the  target  replaced  by  the  free  space  in  a  separate 
calculation. 


As  Example  one,  a  jet  consisting  of  perfectly  conducting  material  with  a 
characteristic  dimension  ^  10  m  is  used  as  the  target.  The  computational 
grid  system  is  shown  in  Fig.  4  with  the  radial  grid  size  Ar  -  In  in  fly 
The  computational  time  increment  At  is  chosen  to  be  1.25ns.  A  dipole  point 
source  with  Gaussian  distribution  in  time  and  a  duration  of  lOAt  is  placed  in  fl^ 
with  the  location  marked  by  "x"  in  Fig.  4.  This  current  source  generates  a 
short  C-W  pulse  of  E-M  fields  with  frequency  *v»  40  MHz  (wavelength  a*  7.5  m). 

The  numerical  results  of  radar  cross  sections  E(6°),  E(16°),  E(35°),  £(62.5°) , 
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and  £(88.5°)  as  functions  of  0  are  plotted  in  Figs.  4,  5,  6,  7,  and  8, 
respectively.  Moreover,  the  radar  back-scattering  cross  section  as  a  function 
of  0(0  £  0  <  90°)  is  plotted  in  Fig.  9;  since  it  is  symmetric  with  respect  to 
9=0,  its  values  in  the  range  of  -90°  <  0  £  0  are  omitted  in  Fig.  9. 

As  a  matter  of  fact,  the  newly  developed  computer  code  in  this  project  is 
used  to  calculate  the  values  of  the  radar  back-scattering  cross  section  in  the 
range  of  negative  9  and  it  is  found  that  the  maximum  deviation  between  its 
corresponding  values  at  either  side  of  0  =  0  is  less  than  3.5%. 

As  Example  two,  a  composite  airfoil  with  the  leading  edge  consisting  of 

-12  -12 

anisotropic  lossy  dielectric,  e  =  22.989x10  farad/m,  e  =  44.210x10  farad/m 
2  x  y 

and  a  =  5x10  mho/m,  ^nd  the  trailing  edge  consisting  of  perfectly  conducting 

material  is  used  as  the  target.  The  characteristic  length  of  the  airfoil  is 

10  m.  The  computational  grid  system  is  shown  in  Fig.  10  and  the  computational 

time  increment  At  is  chosen  to  be  0.15ns.  A  dipole  point  source  with 

Gaussian  distribution  in  time  and  a  duration  of  2.1ns  is  placed  at  (14m, ty) . 

This  current  source  generates  a  short  C-W  pulse  of  E-M  fields  with  frequency 

^  238  MHz  (wavelength  'v  1.26m).  The  numerical  results  of  radar  cross  sections 

£(32.5°),  £(22.5°),  £(12.5°)  and  £(2.5°)  as  functions  of  0  are  plotted  in 

Figs.  11,  12,  13,  and  14  respectively.  Similarly,  the  radar  back-scattering 

cross  section  as  a  function  of  positive  9  is  plotted  in  Fig.  15;  again  due  to 

the  symmetry,  the  radar  back-scattering  cross  sections  for  the  negative  values 

of  0  are  omitted  here. 


TECHNICAL  DISCUSSION 

Obviously,  the  numerical  results  for  the  radar  cross  section  obtained  from 
the  calculation  in  the  previous  section  are  correct  only  in  a  very  approximate 
sense.  The  main  reason  is  that  the  supposedly  scattered  far  field  is 
computed  at  the  location  about  few  wavelengths  away  from  the  target  surface 
(actually  in  the  intermediate  field  region) ,  an  error  in  physical  representation. 
The  simplest  way  to  improve  the  accuracy  is  to  add  many  radial  zones  in  so 
that  the  boundary  will  be  in  the  truly  far  field  region.  Unfortunately, 
this  approach  will  not  really  solve  the  accuracy  problem,  because  the  ratio  of 
the  maximum  dimension  to  the  minimum  dimension  of  the  grid  zone  in  the  far 
field  region  will  be  so  enormous  that  it  will  introduce  large  numerical  errors. 
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Fig.  9  Radar  Back-Scattering  Cross  Section 

To  overcome  this  problem,  a  multi-grid  system  should  be  implemented  in  to 
reduce  the  huge  ratio  of  the  grid  dimensions.  For  this,  an  interpolation 
scheme  must  be  devised  to  transmit  the  information  of  E-M  fields  across  the 
boundary  separating  the  finer  grid  system  from  the  coarser  grid  system. 

There  is  another  computational  efficiency  problem  existing  here.  For  a 
stable  calculation,  the  smallness  of  the  dimension  of  the  triangular  grid  cones 
surrounding  the  origin  of  the  coordinates  makes  At  extremely  small,  and  hence 
it  takes  too  many  At  to  reach  a  p re- determined  time.  To  overcome  this 
difficulty,  one  can  replace  the  cluster  of  small  triangular  grid  zones  by  a 
single  circular  zone  and  modify  the  finite  difference  equations  in  this 
neighborhood  accordingly. 


10  The  shaded  area  represents  the  composite  airfoil 
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xl  -  ( ( float )fabs( theta [ jl] -theta | j] )*x2  + 

(  float)  fobs ( tlietal  j)  -  theta l  j2]  )  *x3)/2,  */ 
xl  -  subt(  theta ,  jl ,  j )  ,■  x5  -  subt.(  theta ,  j  ,  j2 )  ; 
xl  -  <x4  *  x2  +  x5  *  x3  )/2. 0, 
return  (xl)  ,■ 


float  upperr(r,i, j) 
float  r ( ] [ J_MAX J ; 


f  loat  eps_y_i  ( <?pc_y ,  1 
float  eps_y [ ) |J  MAX] , 
int  i,j; 


/*  start  mainline  routine  */ 
int  i, j,k,l,m,n,exchrow 
float  oxchval .innxval ; 


s'S  1 

o  C  t 

'D  — 

**-*  -  X 

ij  t  Ci  ^ 
WsJ  N 

u  « 

-.-4  M  ?3 

^  a? 
X  3  0  0 

O  P  P 
•H  H  W  1 

“5  !C*o  D 
c  x 
«■  w  c  p 
_  'O  ^ 

~  t  CL  0/ 

5  o  o  *a 

P  *4  P  -H 

w  0 

5  C  y  M 


3  S.  'U 
*—  0 
p  »->  o 

-  a,  — • 

3  —  a  **-i 
—  S  H  fl  I  # 

:  •;  C  O  Hi 

-  x  T* 

i.  M  4J  r-H 

*r  o 

J  C  ic 

~  _  S  w  JJ 
3  N 

•  •  —  —•  *H  *H 

g  t  t  <e  «  «  s 

~  “  3  TJ  ^3  P  a 
—  T  T-  t  t  -H  O 

p  c  Copco 

C  -  -  P  P  -H  «P 
<  —  r^n^m 


H  H  H  H  p 
t  O  O  —4  O  p  U) 

P  I  nj  P  I  | 

O  X  O  CO 

p  tj  v  'D  n  n 

t^nO-U*HO»-‘C 
^nOPtOPFt 
3  1  H  I 
0X530S^C 
•H  O  O  O  I  -h 
<t  p  CHiJ  q 

“SiSgJj&S 


s: 

1 

1 

X 

*  « 

3 

« 

C  < 

0  •— 4 

0 

C 

♦-H 

1 

0 

*  « 

rH  C 

« 

-H  -O 

•H  ('r 

1 

T> 

+ 

X 

H 

c 

*  * 

O  *H 

«’ 

U  4p 

-H  P 

3 

0  <6 

c 

a 

< 

cs: 

H 

«  * 

c 

*: 

710 

*4-1  |  3 

E 

P  O 

40 

<a 

_ 

5; 

X 

M 

«  * 

•H  3 

« 

0  F  x 

a 

<4 

t  P 

p 

— 

1 

r* 

0 

*  *  \ 

£ 

*’ 

a  • 

r-4  a  P 

- 

It 

0 

-J  O 

O' 

»M 

«  «;  * 

M  iw 

«■ 

3  *  * 

—4  M  3 

O 

a  * 

P 

3 

•-  O' 

«  « 

« 

/ 

/ 

:  «6  :  0 

1 

S  ^ 

0 

p 

—  0 

0 

t  3 

it 

«  «  w 

'  X 

«r 

si 

a 

0 

•H 

>■ 

— ‘  4-4 

U 

P  0 

p 

«  «  o 

0  P 

« 

X 

C  T> 

a 

•H 

P 

>a 

t  H 

t 

«  «  *H 

•H  P 

* 

O  P  •'  16 

0 

P 

i  tj 

P 

'O  X 

'O 

*  «  -H 

*5  G 

« 

0 

a  G  — 4  O 

Vfl 

<6 

3 

p 

a  • 

p 

«  «  *4-1 

?  x 

p 

0  -H  CLP 

O 

•P 

p 

“ 

iLm  d; 

’O 

0 

0 

«  « 

« 

Q  it  : 

<0 

■H 

C  C-H 

H 

>  'O 

> 

«  ~ 

o 

•  « 

«  ap  « 

E 

H  ■'• 

a 

p 

to-  C"! 

«  Vtol 

0 

It  4-H 

t 

«  «'  X 

c  0 

\ 

\  (t)\ 

O' 

a-— 

a 

3  — 

S  1 

W  0 

a 

«  *. 

:  -H  C 

«• 

0  X 

•H 

p  -h 

0 

p 

«  *H  X 

<M  •» 

«  «'  « 

■  p 

* 

0  p 

« 

p  0 

p 

■A 

1  u 

O  ^  x' 

o  — 

«  *  z 

3  «i-i 

0 

P  -H 

:  > 

T> 

a 

— *  >,-o 

x  a 

pi  a  P 

*  « 

P  i; 

H 

•»3  < 

N 

><V 

'O 

<0 

—  r\ 

P  »-5  c 

•H 

0  -H 

Cu 

«  *  c 

«: 

+ 

:  C*  Q 

J 

C  *H  - 

<6 

p 

— *  x 

i-4  t 

H'O'O  n  T3  CJ 

*  *:  *H 

X 

CL 

>  m  '  . 

3 

5 

0 

4  *  p 

>  ::  S' 

i  :~a 

-  *;  > 

*  O'  c 

J  «'  h  *H  <— 
»  «;  it  <6  — 

!  :  o § 

*  P'S  o  q 

*  :  P  -H  O'  <3 
*:  «0  O'  P  *. 

S:  c'JS  *  P 

*  -H  P  ft 

«  «:  c  x 


POP  \  ^  «J  4—1 

ft  X  -H  ~  *  >,  ft  >  P  O 
•-P  *6  rH  |+J  O  O  > 

U  O  ■'  l/l  Q)  £  4J  (D 

O'  M  ^  >  -v  ai  :  W3  »-H 

S  3%  3  3^1  ~>! 

~  c  «*>  ~  o  x,  p  o=  u* 
OO'O+JX  :  ~  P'U  P 
^4  =  :  <*  3  P  a  'O  =  :  <*  O 
•H  0,  -  a  n  ~  =  C 

.  'w  *M  :|  P  <D  P  'M  — '  O 

f  -M  ^  . — I - 0>P  **-*  | 

P  C  C  O  *0  04*0  C  C4J 
:  -H  -H  t  >  *H  0  «  -H  «t  -H 
C^OOMVhOMOC 
>h  Q.  a  «-t  O'  04  P  Qj  M-H 


\  *  4S  (C  H 

*  >,  re  >  4J  o 

Hp  o  o  > 

-  W  0  X  P  o 


:  p  <u  =  b 

l~  P  TJ  P 


O  P  <0 
5  *H  P 
3  —  O 

0  16  X 

a  p  p 

O  -H 

3  5a 

p  •*.  c 
h  ai  -h 

M  4- 

P  *H  Qi  P 
0  <6  Qj 

O'  I  p  T> 
0  (0 

^&5J! 


int  start , end; 
float  named  [J_MAX] ; 


:h is.  routine  reads  all  levels  on  which  CHEN  wants  to  dump  his  E  and 
[  fields 


int_energy [i] 
hal£_energy [i] 


